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Abstract 

Critical properties of the Ising model on a stacked triangular lattice, with 
antiferromagnetic first and second-neighbor in-plane interactions, are studied 
by extensive histogram Monte Carlo simulations. The results, in conjunc- 
tion with the recently determined phase diagram, strongly suggest that the 
transition from the period-3 ordered state to the paramagnetic phase remains 
in the xy universality class. This conclusion is in contrast with a previous 

suggestion of mean-field tricritical behavior. 
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I. INTRODUCTION 



There is little consensus in recent literature regarding critical phenomena associated 



spin models, noncolinear magnetic ordering gives rise to nontrivial symmetry considera- 
tions and this has lead to contrasting proposals on the nature of the temperature-driven 
phase transition from the period-3 state to the paramagnetic phase: The suggestion by 



put forth by Azaria et alp where first-order, mean- field tri critical or 0(4) criticality can 
occur depending details of the model (also see Refs. [ Even the Ising model on 

this frustrated lattice has a highly degenerate ground statell! and unusual temperature and 
magnetic-field induced phase transitions.^ The conclusion from symmetry arguments is that 
the paramagnetic transition belongs to the standard xy universality class,0 consistent with 
preliminary Monte Carlo result sill More recently, this simple picture has been challenged by 
Heinonen and PetschekEl (hereafter referred to as HP) who made the remarkable proposal 
that this transition is mean-field tricritical, as inspired by their Monte Carlo analysis of 
critical exponents, the structure factor, and the observation that a sufficiently large value 
of third-neighbor in-plane interaction J 3 induces a strong first-order transition to a different 
type of order. (Such an idea has also found support from recent experimental results.0) 
It thus appears that, independently, and for completely different reasons, the suggestion of 
tricriticality associated with the stacked triangular antiferromagnet was made by both HP 
(for the Ising model) and by Azaria et al. (for the Heisenberg model). It is the purpose of 
this work to re-examine the scenario proposed by HP through means of more extensive (and 
more accurate) conventional and histogram Monte Carlo simulations of the Ising model with 
antiferromagnetic first and second-neighbor in-plane interactions (Ji, J2 > 0). 

A number of studies have been made of the effects of further-neighbor interactions on 
the triangular antiferromagnet .EUll At a critical value J2 > 0, the period-3 state (C3) is 
destabilized in favor of magnetic order with a periodicity of two (C2). Of particular interest 



with the simple stacked triangular antiferromag In the 



cases of Heisenberg and xy 
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to the present work are the results of recent Monte Carlo determinations of J 2 — T (temper- 
ature) phase diagrams for IsingEl as well as xy and HeisenbergEl models. The paramagnetic 
transition temperature T/v decreases sharply with increasing J 2 until the critical value, af- 
ter which Tjv increases sharply With J\ = 1, the critical values of J 2 were found to be 
approximately 0.10 for the Ising model and 0.125 for xy and Heisenberg models. For each 
model, the transition line to the C3 state was found to be continuous, within the accuracy 
of these conventional simulations, whereas the C2 transition is always first order. In the 
case of Heisenberg and xy models, a continuous transition to an incommensurate (IC) order 
is observed at larger values of J2 (e.g., J2 ~ 1.1 in the Heisenberg case). 

A simple mean-field argument explains the principal features of these phase diagrams, 
and also reveals a connection with the results of HP. For systems governed by the Hamilto- 
nian 

H = ^^>2 J ijSi ■ s i> (!) 

the Fourier transform of the exchange interaction J(Q) determines the wavevector which 
characterizes the spin modulation of the first ordered state to stabilize as the temperature 
is lowered. (This result is independent of the number of spin components.) In the case of 
ferromagnetic interactions along the c axis (Jo < 0), Q\\ = and the modulation occurs 
entirely in the basal plane. We consider here J(Qj_) with up to third-neighbor in-plane 
interactions includedS Wavevectors which maximize this function give the desired result, 
shown in Fig. 1. With J 3 = 0, the C2 phase is stabilized for | < J 2 < 1 and the IC 
state is realized for J 2 > 1. As found by HP, the C2 state also occurs with J 2 = and a 
ferromagnetic third-neighbor interaction. Their critical value J3 ~ —0.08 can be compared 
with the present mean- field result J 3 = — | ~ —0.111. The scenario put forth by these 
authors calls for interactions between primary C3 and secondary C2 fluctuations to drive 
the transition first order (to the C3 phase) at very small values J3 < 0. Thus J3 — would 
be a tricritical point (not to be confused with the multicritical point where C3, C2 and 
paramagnetic phases meet). In view of the phase diagram Fig. 1 and the results of Ref. [ 
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18| , such behavior should also be revealed by considering the effects of small J 2 (as done 
in the present study). Finally, note that the dramatic depression in T/v at the multicritical 
point observed in the Monte Carlo results of Ref. [ [18] is seen from Fig. 1 to likely be a 
result of large critical fluctuations since the IC and C2 phases are degenerate at this value 
of J 2 . The relatively weak dependence of T N on J 3 < observed by HP is consistent with 
this explanation. 

Although the conventional Monte Carlo simulations of Ref. [ |18| suggest that the tran- 
sition with < J2 < 0.1 remains continuous, the implementation of more sensitive finite- 
size scaling techniques may be required to detect a very weak first order transition, as 
might be expected near a tricritical point. The Ferrenberg-Swendsen histogram Monte 
Carlo methodic has proven useful for this purpose, especially when used to determine the 
limiting value of the (internal) energy cumulant@ 

U(T) = 1 - - < E A > I < E 2 > 2 . (2) 
3 

This quantity has a minimum at a phase transition, with the property U(Tn) — > U* — | in 
the infinite-lattice limit. The histogram method allows for the possibility of precise determi- 
nation of extrema exhibited by other thermodynamic functions at T/y. The scaling behavior 
with system size of these quantities can yield accurate estimates of critical exponents in the 
case of a continuous transition, or reveal simple volume dependence if the transition is first 
order. The utility of this type of Monte Carlo method for the present purposes is very nicely 
described by Reimers et a/.@ Recently, Bunker et a/.@ made such histogram analysis of the 
present model in the case of nearest-neighbor in-plane interactions only. They find critical 
exponents consistent with xy universality, in contrast with HP. Our work corroborates and 
extends their results. 



Guided by the results of Ref. [ [Te[] , we test the proposals of HP by performing Monte Carlo 
simulations at three values of J 2 : 0, 0.08, and 0.25. Conventional Monte Carlo analysis at 
J 2 = reveals that the critical region is rather narrow in temperature, thus providing 
an expalanation for the erroneous exponent estimates given by HP. Extensive histogram 



simulations at J 2 = 0.08 provide convincing evidence that the transition remains continuous 
and of xy universality. The results at J 2 = 0.25 serve as an example of scaling behavior in 
the case of a first order transition to the C2 state. 



II. J 2 = 

Conventional Monte Carlo simulations were performed for the case of nearest-neighbor 
in-plane interactions only. Runs of 2xl0 4 to 10 5 Monte Carlo steps (MCS) per spin were 
made with the initial 4xl0 3 to 2xl0 4 MCS discarded for thermalization. In the case of 
smaller runs, quantites were averaged over 4 independent simulations using random initial 
spin configurations. Periodic boundary conditions on lattices LxLxL with L=12-30 were 
used. Finite-size scaling of only the critical exponent j3 was considered, using the C3 order 
parameter M ~ t@ [t = (Tv — T)/Tv)] defined in terms of a Fourier component as in Ref. [ 
28"f. The extrapolation technique of LandauH was used on results close to the transition 



temperature, known from HP to be near 2.9Ji, to estimate their values for L — > 00. The 
results are presented in Fig. 2, where L=12 data were excluded in the fit for the highest four 
temperatures. Tv was then adjusted to yield the best linear fit of the extrapolated data for 
/n(Moo) vs. ln{t) plots. 

Initial analysis performed on data for 12 < L < 24 and 2.0 < T < 2.7 yielded values 
Tjv — 2.82 and (3 ~ 0.19, close to those of HP. However, a more detailed study using data 
at larger L and higher T revealed a different set of results. It was difficult to achieve good 
linear fits of all the data between T=2.00 and T=2.90 so that some of the lower-temperature 
points were excluded. The outcome of this procedure was instructive. Using only the six 
highest temperatures (T=2.60-2.90) gave the estimates Tv — 2.92 and (3 ~ 0.30. Fig. 3 
shows results with only the four highest temperatures included, which yield the estimates 
Tv — 2.93(1) and (3 ~ 0.33(3). The latter values are dramatically different from those of HP 
but are consistent with more accurate histogram results (see below and Ref. [ ^6|) and the 
expected xy universality. The above analysis strongly suggests that this model exhibits an 
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unusually narrow critical region of temperature. This may be a consequence of the proximity 
of the C2-order instability and the coupling of these fluctuations to the primary C3 order, 
as suggested to be important for this model by HP. 

Such a system is ideally suited to exhibit the power of the histogram method, where 
exponents can be extracted from simulations performed very close to the transition tem- 
perature. Since these simulations were done mainly as a test for later runs at J2 = 0.08, 
and to corroborate the more detailed study of Ref . [ EJf| , only a single histogram was made 
for each lattice size L=12,15, 18,21,24,30 at the temperature T=2.93 with 1.2xl0 6 MCS per 
spin, where 2xl0 5 MCS were discarded for thermalization. More details are presented in the 
next section for the case of J2 = 0.08 and, in the interest of brevity, we present no data here 
but state only the results of our analysis. Finite-size scaling of the maxima (or minima) for 
the specific heat (C), susceptibility (x), energy cumulant (U), and logarithmic derivative of 
the order paramete 

V(T) = < ME > I < M > - < E >, (3) 

yield estimates for a/u, 7/z/, U*, and 1/z/, respectively. In addition, M{T N ) scales as L - ^, 
where T/v can be adjusted to give the best linear fit. The results a/u < 0.2(3), (3/v = 0.50(3) 
(using T N = 2.928), 7/z/ = 2.03(6), 1/z/ = 1.46(5), and U* ~ 0.666667(5) were obtained. 
(Errors are difficult to assign; those given here are estimated from the robustness of the linear 
fitting. Evaluation of results from many simulations are necessary to obtain more reliable 
estimates of statistical errors). These values compare favorably with renormalizat ion-group 
results a/u = -0.018, f3/i> = 0.519, 7/z/ = 1.96, l/v = 1.49. The quantity (U* - U) scales 
as L~ m , where in principle m = a/u. However, only very weak maxima in C and minima 
in U were observed,!^ consistent with a very small value of a. As noted in Ref. [EM], in such 



cases (U* — U) may exhibit simple volume dependence. The observed exponent m c± 2.75 is 
consistent with this expectation. 
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III. J 2 = 0.08 



Finite-scaling analysis of histogram Monte Carlo simulations for the case of J 2 = 0.08 are 
presented here. A relatively large value of second-neighbor coupling was desired to enhance 
the possibility of observing first-order effects, but not too close to the multicritical point 
at J 2 ~ 0.1 where larger fluctuations might be expected. However, since J 2 is quite small 
relative to the other energy scale in the model (Ji = 1), longer runs are necessary to fully 
realize its effects (and achieve the same accuracy as with J 2 = 0). With these considerations, 
it can be expected that the extraction of critical exponents will be less reliable than for the 
case of J 2 = 0. 

Simulations were performed using the same parameters as described above for the case 
J 2 = 0, but at two or three different temperatures near the transition estimated from 



the results of Ref.[ [18[ to be T/v — 2.2. This ensured that the histogram was made at a 
temperature close enough to the extrema of the thermodynamic quantity of interest for the 
extraction of meaningful results. As a general guide to reliable data, the temperature of the 
extrema should lie within the range determined by half of the maximum of the histogram. 
If it does not, a new histogram should be generated at a temperature closer to the extrema 
of the desired function. 

Logarithmic scaling plots were made of the various thermodynamic functions as described 
in the previous section. Convincing evidence that the transition remains continuous is found 
in results for the energy cumulant, where the value U*=0. 666655(20) was extracted. Fig. 4 
displays results for (U* — U), with the assumtion U* = |; the good linear fit further supports 
this conclusion. The specific heat (as well as U) again exhibited only very weak temperature 
maxima, with a small dependence on L, as displayed in Fig. 5. It can be concluded from 
these data that ajv is very small. Scaling of the other functions produced the estimates 
P/v = 0.50(3) (using T N = 2.197), yjv = 2.07(6), \jv = 1.40(5), and m = 3.10(10), 
obtained from all of the data. Using data from only the three largest lattice sizes yielded 
slightly different results for y/v = 2.00(6), 1/v = 1.42(5), and m = 3.04(10). In an effort 
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to demonstrate further that our data are consistent with xy universality, the scaling plots 
shown in Figs. 6-9 were made using the renormalizat ion-group exponents. The results are 
convincing on this point. We note that no visual discrimination could be made between 
these plots and those made with our extracted exponents. 



IV. J 2 = 0.25 

Finally, a brief summary of histogram results for the case of J 2 = 0.25 where a rela- 
tively strong first-order transition to the C2 phase is expected.lii The main purpose of these 
simulations is to illustrate that the energy cumulant is sensitive to the order of the transi- 
tion for the present model. Since the ordered state has a periodicity of two, even values of 
the lattice size L were used. Simulations were perfomed with L=12,18,24,30 at a range of 
temperatures guided by considerations as outlined in the previous section, near the known 
transition point, T/v — 2.66. Very large fluctuations were found in the results for the x an d 
V, likely a consequence of metastability effects. Finite-size scaling of these data was not 
fruitful. Fig. 10 displays scaling results for the minima in U, which were used to obtain 
the saturation value U*=0. 6595(3). This result, along with a reasonably good linear fit of 
the specific heat maxima vs. L 3 , confirms the strong first-order nature of this transition and 
further strengthens the conclusions of the previous section. 



V. SUMMARY AND CONCLUSIONS 



The results presented in this work and in Refs.[ pl| and |26|] strongly suggest that the 
phase transition associated with the period-3 state of the next-nearest-neighbor Ising an- 
tiferromagnet on a stacked triangular lattice exhibits xy universality. In contrast with the 
proposal by Heinonen and Petschek, no mean-field tricritical behavior is observed. The idea 
by these authors that the coupling of fluctuations to the nearby period-2 order should be 
important, however, appears to be responsible for a shrinking of the true critical region 
and offers an explanation for their erroneous estimates of critical exponents. The histogram 
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Monte Carlo method was found to be ideally suited for the exposition of the real critical 
behavior in this system. It remains to be tested if these conclusions are relevant in the 
case of continuous spin models on this frustrated lattice where there are diverse proposals 
regarding the criticality. 
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FIGURES 

FIG. 1. Phase diagram determined by maximizing J(Qj_) with up to third-neighbor in-plane 
interactions (with antiferromagnetic first-neighbor coupling J\ = 1) where C2 and C3 represent 
commensurate phases of periodicity two and three, respectively, and IC denotes the incommensu- 
rate phase. Solid and broken lines indicate first and second-order transitions, respectively. 

FIG. 2. Finite-size scaling of the order parameter at J2 = and selected temperatures 
T=2.00-2.90 for lattice sizes L=12-30. L=12 values were excluded from the fit for the four highest 
temperatures. 

FIG. 3. Ln-ln plot of the order-parameter values for J2 = extrapolated to L — > 00 from the 
results of Fig. 2. Only the four highest temperatures were used in the linear fit, with the results 
T N ~ 2.93(1) and (3 ~ 0.33(3). 

FIG. 4. Scaling of the energy-cumulant minima for J 2 = 0.08 with U* set to §. The resulting 
estimate for m is given by the slope. 

FIG. 5. Temperature behavior of the specific heat showing the weak maxima and small lat- 
tice-size dependence for L=12 (lower curve) to L=30 (upper curve). 

FIG. 6. Scaling behavior of the specific heat maxima assuming the renormalization-group value 
of a/v for the xy universality class. 

FIG. 7. Scaling behavior of the order parameter at Tjv = 2.195 assuming the renormaliza- 
tion-group value of (3/v for the xy universality class. 

FIG. 8. Scaling behavior of the susceptibility maxima assuming the renormalization-group 
value of 7/1/ for the xy universality class. 

FIG. 9. Scaling behavior of the maxima in the logarithmic derivative of the order parameter 
(3) assuming the renormalization-group value oil/v for the xy universality class. 
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FIG. 10. Scaling of the energy-cumulant minima with L for the case J2 = 0.25. The extrapo- 
lated estimate for U* is 0.6595(3). 
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